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ABSTRACT 



The discovery of bright gamma-ray emission coincident with supernova rem- 
nant (SNR) W51C is reported using the Large Area Telescope (LAT) on board 
the Fermi Gamma-ray Space Telescope. W51C is a middle-aged remnant (~ 10^ 
yr) with intense radio synchrotron emission in its shell and known to be inter- 
acting with a molecular cloud. The gamma-ray emission is spatially extended, 
broadly consistent with the radio and X-ray extent of SNR W51C. The energy 
spectrum in the 0.2-50 GeV band exhibits steepening toward high energies. The 
luminosity is greater than 1 x 10^^ erg given the distance constraint oi D > 5.5 
kpc, which makes this object one of the most luminous gamma-ray sources in our 
Galaxy. The observed gamma-rays can be explained reasonably by a combination 
of efficient acceleration of nuclear cosmic rays at supernova shocks and shock- 
cloud interactions. The decay of neutral vr-mesons produced in hadronic colli- 
sions provides a plausible explanation for the gamma-ray emission. The product 
of the average gas density and the total energy content of the accelerated protons 
amounts to n^Wp ~ 5 x 10^^ (-D/6 kpc)^ erg cm~^. Electron density constraints 
from the radio and X-ray bands render it difficult to explain the LAT signal as 
due to inverse Gompton scattering. The Fermi LAT source coincident with SNR 
W51G sheds new light on the origin of Galactic cosmic rays. 

Subject headings: acceleration of particles — ISM: individual(W51G) — radiation 
mechanisms: non-thermal 
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Introduction 



The o rigin of cosmic r ays remains unsolved. The idea that supernovae produce cosmic 



rays (e.g., iHayakawal Il956l ) has been developed both in theoretical and observational as- 
pects, so that galactic cosmic rays are widely considered to be accelerated in the expanding 
shocks of supernova remna nts (SNRs) through diffusive shock acceleration process (see e.g., 
Blandford fc Eichled 119871 ). This conjecture has been strengthened by recent observations 
of young SNRs in synchrotron X-rays and very-high-energy (VHE) gamma-rays (see e.g.. 



Reynolds! l2008l : lAharonian et al.ll2008l ). High efficiency for converting the kinetic energy 
released by supernova explosions into the energy of relativistic protons and nuclei is a key 
requirement to explain the galactic cosmic rays; it is yet to be confirmed observationally. 

Enhanced vr^-decay emission expected from SNRs that are interacting with molecular 
clouds could provide direct evidence of the nuclear component of cosmic rays being accel- 
erated at supernova shocks (lAharonian et al.lll994j ). To identify the vr^-decay gamma-rays 
as evidence of the nuclear cosmic rays, observations in the GeV domain are crucial because 
the TT^-decay spectrum has a characteristic bump around 70 MeV. Previous measurements 
of GeV gamma-rays with EGRET onboard the C ompton Gamrna-Ray Observatory found 
some gamma-ray sources near radio-bright SNRs (lEsposito et al.lll996l ). However, the pos- 
sible origins of the EGRET sources, namely SNR shell contributions or pulsar associations, 
remained unclear, mainly due to poor localization. 

The advent of the Large Area Telescope (LAT) onboard the Fermi Gamma-ray Space 
Telescope provides a new opportunity to study the gamma-ray emission from SNRs at 
GeV energies. An initial source list based on the first three months of Fermi observa- 



tions (lAbdo et al.ll2009d ) includes OFGL J1923. 0-1-1411, which is spatially coincident with 
SNR W51C. Even with the three months data, the detection was at ~ 23a level. There are 



no EGRET counterpart (s) to the LAT source (IHartman et al.l 119991 ): this is the first SNR 
discovered by Fermi as confirmed in this paper. 

W51C (G49.2— .7) is a radio-bri ght SNR at a distance of D ~ 6 kpc with an estimated 
age of ~ 3 X 10^ yrs (iKoo et al.lll995l ). It has an elliptical sh ape with an extent of 50' x 38 ^ 
The radio continuum map exhibits thick shell-like structures (ISubrahmanyan &: Gosslll995l ). 
The bul k of the X-ray em ission comes from thermal plasma with a temperature of kT ~ 
0.3 keV (IKoo et al.ll2005l ). A molecular cloud-shock interaction is known , as evidenced by 
observations of shocked atomic and molecular gases (IKoo fc Moonlll997al Jbl). Most recently, 
the HESS collaboration has announced the detection of extended VHE gamma-ray emission 
coincident with W51C (IFiasson et al.ll2009l ). Also, the Mil agro collaboration h as reported a 
possible excess of multi-TeV gamma-rays in this direction (lAbdo et al.ll2009af ). 
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In this paper we report the analysis results for the LAT source coincident with SNR 
W51C, using data accumulated over the first year of Permits operation. The Fermi observa- 
tions of SNR W51C permit a refined study of cosmic-ray acceleration. Specifically, the LAT 
data suggest that 7r°-decay emission is the dominant contribution to the gamma-ray signal. 



Observation and Data Reduction 



The Fermi Gamma-ray Space Telescope was launched on 2008 June 11 by a Delta II 
Heavy launch vehicle. The LAT onboard Fermi is a pair-conversion gamma-ray detector 
capable of measuring gamma-rays in a very wide range of energy from 20 MeV up to 300 
GeV. The LAT tracks the electron and positron resulting from pair conversion of an incident 
gamma-ray in thin high-Z foils, and measures the energy deposition due to the subsequent 
electromagnetic shower that develops in the calorimeter. The effective area is ~ 8000 cm^ 
above 1 GeV (on-axis) and the per-photon 68% containment radius is ~ 0?8 at 1 GeV. 
The point-spread function (PSF) depends largely on photon energy and improves at higher 
energies. The tracker of the LAT is divided into two regions, front and back. The front region 
(first 12 planes) has thin converters to optimize the PSF while the back region (4 planes after 
the front section) has thicker converters to enlarge the effective area. The angular resolution 
for the back events is approximately twice as broa d than that for the front events. The 
details of the LAT and da ta processiiig are gi ven in lAtwood et al.l ( 120091 ) . and the on-orbit 
calibration is described in lAbdo et al.l (l2009bl ). 



The gamma-ray data acquire d from 2008 August 5 to 2009 July 14 are analyzed. The 
diffuse class events as defined in lAtwood et al.l (120091 ) are chosen for gamma-ray analysis. 
A cut on earth zenith angles greater than 105° is applied to reduce the residual signal from 
earth albedo gamraa-ray s. The instrument response functions (IRFs) called "Pass 6 V3" are 
used jRando et al.ll2009h . 



3. Analysis and Results 

The maximum likelihood technique is employed for spectral and spatial parameter es- 
timation using gtlike, which is publicly available as part of Fermi Science Toolfl The like- 
lihood is the product of the probability of observing the gamma-ray counts of each spatial 



^Software and documentation of the Fermi Science Tools are distributed by the Fermi Science Support 
Center at |http: / / fermi.gsfc.nasa.gov/ssc, 
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and energy bin given the emission model, an d parameter values a re estimated by maximizing 
the likelihood of the data given the model (IMattox et al.lll996l ). The gamma-ray emission 
model includes individual sources at fixed coordinates, galactic diffuse emission (resulting 
from cosmic-ray interactions with interstellar gas and photons), and an isotropic component 
(extragalactic and residual background). The so-called "mapcube" file of glLiem_v02.fit is 
used for modeling the galactic diffuse emission, together with the corresponding tabulated 
model for the isotropic diffuse emission. Other v ersions of the galactic diffuse models, gener- 
ated by the GALPROP code ( IStrong et al.ll2004j ). are also utilized to assess systematic error. 
The maximum likelihood analysis is performed inside a square region of 12° x 12° centered 
on W51C with a pixel size of 0?1, unless otherwise mentioned. Background point sources 
detected in six months data are included in the likelihood analysis with free normalization 
and power-law index, though none of them affect the results in this paper. 



3.1. Spatial Distribution 

In Figure [U the maps of photon counts in the 2-10 GeV band in the vicinity of SNR 
W51C are shown; the right panel is a close-up view of the left panel. Gamma-ray events 
that converted in the front section of the tracker are selected. A bright gamma-ray source 
is enclosed by the outer boundary of W51C. The average surface brightness inside the SNR 
boundary is about 2 and 5 times larger than neighboring regions on the galactic plane in the 
0.5-2 GeV and 2-10 GeV bands, respectively. The gamma-ray distribution appears to be 
somewhat clumpy, suggesting the presence of sub-structures. Due to the limited statistics, 
we defer the investigation of possible sub-structures to a future publication. 

The spatial distribution of the gamma-ray source is significantly extended compared to 
a simulated point source that has the same spectrum. A one dimensional profile in the right 
ascension direction of the counts map is compared with that expected for a point source 
in Figure [2l Though the width of the PSF above 5 GeV is known to deviate from the 
Monte Carlo simulations, it has negligible effects on the simulated point source and other 
results in this paper. Assuming a two-dimensional gaussian distribution fixed at the W51C 
centroid (a, (5) = (290?818, 14? 145), we calculate the extension parameter of the source as 
(Text = 0?22 ± 0?02 by varying Uext to find the best match with the data. Note however that 
a different assumption on the spatial distribution results in different quantification of the 
spatial extension. 

The extended nature of this LAT source cannot be ascribed to a superposition of two 
point-like sources. To quantify this, we define a grid of 60 positions inside the remnant 
as trial point source positions, and calculate maximum likelihood values {L2ps) for each 
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possible combination by running gtlike for a box of 8° x 8°. On the other hand, assuming 
that the gamma-ray source has uniform surface brightness inside the SNR boundary, we 
obtain another maximum hkehhood value (Lum)- The resulting values of the likelihood test 
statistics, — 2 ln(L2ps/-^vuni) > 32, demonstrate that the uniform emission gives a much better 
result. 

The extension of the source precludes magnetospheric radiation from a pulsar as a 
dominant gamma-ray source. However, a small fraction (~ 10%) of the gamma-ray flux 
could be contributed by a pulsar. Our pulse search results in non-detection, which places a 
5(7 upper limit on the p ulsed gamma-rays as ~ 6 x 10~^ photons cm~^ s~^ above 100 MeV 



see 



Abdo et al.ll2009dl ) 



3.2. Spectrum and Its Uncertainties 



The gamma-ray spectrum of the W51C source is shown in Figure O It is obtained by 
performing likelihood analysis for each energy bin with gtlike. The lower energy bound is set 
at 0.2 GeV, below which the systematic uncertainties become too large due to the background 
diffuse emission (see below). The surface brightness of the W51C source is assumed to be 
uniform inside the SNR boundary. The normalization of the galactic diffuse emission model 
is left free at each energy bin. As evident from Figure [31 the gamma-ray spectrum cannot 
be described by a simple power law and steepens above a few GeV. A likelihood-ratio test 
indicates that a power-law hypothesis has a chance probability of 5 x 10~^ for obtaining 
the spectral data. The physical interpretation of the spectral energy distribution will be 
discussed in §11 

The gamma-ray luminosity can be estimated as ~ 1 x 10^^ {D/6 kpc)^ erg s~^ in the LAT 
bandpass, making this object one of the most luminous gamma-ray sources in o ur Galaxy. 
Note that the distance to the remnant is well constrained. X-ray absorption (IKoo et al. 



19951 ) indicates that W51C is situated behind the W51 molecular cloud complex (a so-called 
68 km s~^ cloud in particular) w hose distance is de termined as 7.0 ± 1.5 kpc by the proper 
motion of W51 Main H2O mas er (IGenzel et al. I ll98lh . On the other hand, the angular extent 
and the X-ray measurements ( Koo et al. 2002h requi re D < 10 kpc. We adopt D = 6 kpc 
following previous pubhcations (e.g.. lKoo et al.l 120051 ). 



It should be noted that in addition to the systematic errors commonly assigned to LAT 
spectral dataE, the accuracy of the flux estimated for each energy bin of the W51C source 



^For the IRF used here (P6.V3_DIFFUSE), systematic error as a function of a; = log(£;/MeV) is 10% at 
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is limited by the following errors. First, a possible effect of the uncertainty of underlying 
galactic diffuse emission is considered. The observed gamma-ray intensity of nearby source- 
free regions on the galactic plane is compared with the intensity expected from the galactic 
diffuse models. The difference, namely the local departure from the best-fit diffuse model, 
is found to be < 6%. By changing the normalization of the galactic diffuse model artificially 
by ±6%, we estimate the possible systematic error to be 40% (0.2-0.4 GeV), 22% (0.4-0.8 
GeV) and < 10% (> 0.8 GeV). 

The fact that we do not know the true gamma-ray morphology of the W51C source 
introduces another error in our flux estimation. Since the gamma-ray source is point-like 
below 1 GeV given the PSF, this uncertainty matters only above 1 GeV. We adopt a uni- 
form surface brightness inside the SNR boundary (a flat elliptical template) as the spatial 
distribution of the source gamma-rays. Different spatial distributions such as a fiat elliptical 
template reduced in size (scaled by 0.5) are tested to estimate the systematic error. Our 
conservative estimate is < 20% in 1-6 GeV and ~ 30% above 6 GeV as the systematic 
uncertainty attributable to the unknown shape of the source. 



The extended gamma-ray emission positionally coincident with SNR W51C has been 
studied using the Fermi LAT. The gamma-ray spectrum presented in Figure [3] is not fit 
by a simple power law, exhibiting a remarkable steepening. Here we discuss the origin of 
the extended emission and the underlying particle spectra that give rise to the observed 
spectrum of photons. 

The expanding shock waves driven by the supernova explosion are expected to be the 
sites of the acceleration of multi-GeV particles. To phenomenologically interpret the spec- 
tral curvature in the LAT-l-TeV bands, a broken power-law is adopted for the momentum 
distribution of the radiating electrons/protons: 



where Po = 1 GeV c~^. For simplicity, the indices and the break momentum are assumed to be 
identical for both accelerated protons and electrons. As we argue below, the break may reflect 
the character of magnetohydrodynamic turbulence. To account for the radio synchrotron 



4. 



Discussion 




(1) 



X = 2, 5% at X = 2.75, and 20% a,t x — A with a linear interpolation between them. 
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index a ~ 0.26 (iMoon fc Kodll994j ). we adopt s = 1.5, though a steeper index (say, s = 
1.7) could be reconciled with the radio observations within the uncertainty. The energetic 
particles are assumed to be uniformly distributed in the volume of ^ = (47r/3)-R||g with an 
effective radius R^s = 30 pc. The age and radius imply a shock velocity of Vgh ~ 400 km s~^ 
and Est^/rio ~ 1.6 x 10^^ erg cm^, where £^sn and no represent the explosion kinetic energy 
and the interstellar density into which the main blast wave is propagating, respectively. The 
X-ray data suggest uq ~ 0.3 cm~^. The radio images in dicate that rad i o-emit ting electrons 
are smoothly distributed in a thick shell. The model of iKoo fc MoonI (jl997al ) suggests the 
presence of a molecular cloud of ~ 1 x 10^ Mq engulfed by the blast wave. The engulfed 
cloud can act as target material for relativistic particles. The total (atomic and molecular) 
hydrogen mass contained in the volume is denoted by Mh = n-nrripV. Note that ^ uq 
can be expected. 

Given the interaction with a molecular cloud, we first attribute the observed gamma- 
rays to the decay of 7r° mesons produced in inelastic collisions between accelerated protons 
(and nucl ei) and target gas ( Fig. [3]). The gamma-ray spectrum of 7r°-decay is calculated 
based on iKamae et al.l (120061 ) using a scaling factor of 1.85 for helium and heavy nuclei 
( lMorill2009l ). Note that the scaling factor assumes the local interstellar abundance for target 
material and the observed cosmic-ray composition. Contributions from bremsstrahlung and 
Inverse- Compton (IC) scattering by accelerated electrons are also show n in Fig. [3l E l ectron - 
ion and electron-electron bremsstrahlung spectra are computed as in iBaring et al.l (Il999l ). 
The interstellar radiation field for IC (see Table [T]) is comprised of two diluted blackbody 
components (infrared and optical) and the CMB. The infrared and optical co mponents are 



adjusted to reproduce the interstellar radiation field in the GALPROP code (jPorter et al. 



20081 ). Cooling effects due to ionization and synchrotron (or IC) losses, which introduce 
cooling breaks in particle spectra in addition to pbr, are taken into account assuming constant 
particle injection over a period Tq ~ 3 x lO'^ yr. The synchrotron cooling becomes important 
in the TeV band for leptonic models. 

Figure H] (a) shows the radio+gamma-ray spectrum together with the radiation model 
that uses the parameters in Tabled! We adopt here Mh = 2.8 x lO'' Mq (tt-h = 10 cm~^), 
which is somewhat larger than the value quoted above. The total energy of the high-energy 
protons amounts to Wp = 5.2 x 10^° erg, which is inversely proportional to Mh, but insensitive 
to other parameters. The large luminosity of the LAT source can be explained naturally by 
a large number of accelerated protons and dense environments as expected for the case of 
W51C. 



The break feature in the proton spectrum, which is introduced on phenomenological 
grounds, is not expected in a shock acceleration zone if acceleration proceeds close to the 
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Bohm limit (e.g.. iBaring et al.lll999l ). Protons can be accelerated up to ~ 200 TeV via diffu- 
sive shock acceleration with a power-law (or even slightly concave) momentum distribution. 
A possible explanation for the falling proton spectrum above pbr = 10-30 GeVc~^ (which 
depends on the choice of other parameters) could be the effects of damping of magnetohy- 
drody namic turbulence due to ion-neutral collisions. According to iPtuskin fc Zirakashvili 
(120031 ). the maximum attainable momentum in the presence of wave dissipation by ion- 
neutral colhsions is estimated as Pmax(^o) ~ 30 (n/1 cm"'^)^^/^ GeVc^^, using the SNR 
parameters described above. Here n is the ambient neutral gas density. The shock pre- 
cursor cannot confine accelerated protons with p > Pmax(To) even if they were accelerated 
earlier. Therefore, energy-d ependent escape of accelerated particles from the remnant (e.g., 
Aharonian &: Atoyanlll996l ) may account for the steepening of the proton distribution. In 
any case, a break is needed to fit the observed SED. 

In contrast, leptonic models for the GeV gamma-ray production face difficulties (see 
Figure Hb-c and Tabled]). While the chosen model parameters are not unique in their ability 
to fit the broadband spectrum, they are representative; no reasonable choices for them can 
easily account for the radio+gamma-ray data. Leptonic scenarios require Oe/ap far in excess 
of cosmic-ray abundance ratios. It is difficult to reconcile the bremsstrahlung-dominated 
model with the observed radio synchrotron spectrum (Fig. Ht). Also, the IC-dominated 
model (Fig. Ht) requires an unrealistically large energy content in radiating electrons. We — 
1 X 10^^ erg, a low magnetic field of 5 ^ 2 /iG (to suppress the radio synchrotron emission), 
and a low density of nn < 0.1 cm~^ (to reduce the bremsstrahlung component). Such a low 
density is problematic because even X-ray-emitting gas has a density of ~ 1 cm~^ (IKoo et al. 



20021 ). It is also difficult to account for the LAT signal by an IC process in a possible pulsar 
wind nebula seen inside W51C, CXO J192318.5+143035 JKoo et al.ll2005l ). which has only 
a few parsec extent and much weaker radio emission compared with the SNR shell. 

The most plausible explanation for the LAT source is therefore offered by vr^-decay 
gamma-rays in a dense environment, spawned by efficient acceleration of protons and nuclei 
taking place or having taken place at the shocked shell of SNR W51C. The discovery of the 
GeV-scale gamma-rays presented in this paper raises the intriguing possibility that Fermi has 
discerned accelerated ions in a middle-aged remnant that is interacting with dense clouds. 
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Fig. 1. — (Left) Fermi LAT counts map in 2-10 GeV around SNR W51C in units of counts 
per deg^. Front-converted events are selected. The counts map is smoothed by a Gaussian 
kernel of cr = 0?12. The green dashed line represents the galactic plane. (Right) Close 
up view around SNR W51C. The simulated point source image (smoothed by the same 
gaussian) is shown in the inset. The outer boundary of W51C is indicate d by a white ellipse. 
Superposed is the R OSAT X-ray rnap (co ntours) from iKoo et al.l (119951 ). The region where 
shocked CO clumps (IKoo fc Moonlll997bl ) were found is represented by a dashed magenta 
ellipse. A diamond near the SNR centroid indicates CXO J192 318. 5-1-143035 (see the text) . 
The positions of five HII regions are indicated by green crosses (jCarpenter fc Sanderslll998l ). 
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Fig. 2. — One dimensional profile of the observed 2-10 GeV gamma-rays (data points) along 
right ascension. The counts map is integrated over the direction of declination with a width 
of AS = 1?6 centered on W51C. The dashed curve shows the profile of the galactic-|- isotropic 
diffuse model. The histogram represents the sum of a simulated point source and the diffuse 
model. 

Table 1. Parameters of Multiwavelength Models 



Parameters Energetics 



Model cie/ap As pbr B nn Wp We 

(GeV c-i) (/xG) (cm-3) (lO^o erg) (10™ erg) 



(a) TT^'-decay 


0.02 


1.4 


15 


40 


10 


5.2 


0.13 


(b) Bremsstrahlung 


1.0 


1.4 


5 


15 


10 


0.54 


0.87 


(c) Inverse Compton 


1.0 


2.3 


20 


2 


0.1 


8.4 


11 




Note. — Seed photons for IC include the CMB {kTcuB = 2.3 x 10"^ eV, f/cMB = 
0.26 eV cm-3), infrared (/cTir = 3 x 10"^ cV, Um = 0.90 eV cm^^), and optical {kTopt = 0.25 
eV, [/opt = 0.84 eV cm~^). The total energy content of radiating particles, We^p, is calculated 
for J? > 10 McVc^i. 
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Fig. 3. — SED of the SNR W51C region measured with the Fermi LAT (red points) together 
with phenomenological modehng. Systematic errors (see §3.2p are indicated by black bars. 
The model consists of vr'^-decay (long-dashed curve), bremsstrahlung (dashed curve), and IC 
scattering (short-dashed curve). The integrated flux reported by HESS is converted to the 
differential flux at 1 TeV assuming photon index F = 2.5 ± 1.0 (black point). 
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Fig. 4. — Three diff erent scenarios for t he multiwavelength modehng (see Table [T]). The 
radio emission (from iMoon &: Kod Il994l ) is explained by synchrotron radiation, while the 
gamma-ray emission is modeled by different combinations of 7r°-decay (long-dashed curve), 
bremsstrahlung (dashed curve), and IC scattering (dotted curve). The sum of the three 
component is shown solid curve. 



